using Printf

function gen_qrv(mdp, ndx, pro, lig; withLig, qrv, radType, radLJ0)
    qrv = open(qrv, "w")
    write(qrv, "Protein Ligand\n")

    ndxPro = ndx.groups[pro].indexes
    if withLig == true
        ndxLig = ndx.groups[lig].indexes
    end
    
    mdp_file = readlines(mdp)

    # 记录力场参数
    locater = findfirst(x -> occursin(r"\s*ffparams:", x), mdp_file)
    Atyp = parse(Int32, match(r"\s*atnr=(\d+)", mdp_file[locater + 1]).captures[1])
    write(qrv, string(Atyp) * "\n")
    ntypes = parse(Int32, match(r"\s*ntypes=(\d+)", mdp_file[locater + 2]).captures[1])
    locater += 3
    sigma = zeros(Float64, Atyp)
    epsilon = zeros(Float64, Atyp)
    Rad = zeros(Float64, Atyp)
    for i in 1 : Atyp
        write(qrv, @sprintf("%6d", i - 1))
        # 获取每个原子的C6和C12
        for j in 1 : Atyp
            C6 = match(r".*c6\s*=\s*(.*),.*", mdp_file[locater]).captures[1]
            C12 = match(r".*c12\s*=\s*(.*)", mdp_file[locater]).captures[1]
            write(qrv, @sprintf(" %s %s", C6, C12))
            # 计算每种原子类型的σ, ε, 半径
            if j == i
                sigma[i] = 0
                epsilon[i] = 0
                Rad[i] = radLJ0
                C6 = parse(Float64, C6)
                C12 = parse(Float64, C12)
                if C6 * C12 != 0
                    sigma[i] = 10 * (C12 / C6) ^ (1 / 6) # 转换单位为A
                    epsilon[i] = C6 ^ 2 / (4 * C12)
                    Rad[i] = 0.5 * sigma[i]          # sigma为直径
                end
            end
            locater += 1
        end
        write(qrv, "\n")
    end

    Ntyp = parse(Int32, Utils.get_mdparam(mdp, "#molblock"))
    reg = r"\s*moltype\s*=\s*(.+)\s*"
    Imol = parse.(Int32, getindex.(split.(getindex.(match.(reg, filter(x -> occursin(reg, x), mdp_file)), 1)), 1))
    reg = r"\s*#molecules\s*=\s*(.+)\s*"
    Nmol = parse.(Int32, getindex.(match.(reg, filter(x -> occursin(reg, x), mdp_file)), 1))

    reg = r"atom \((.+)\):"
    maxNatm = max(parse.(Int32, getindex.(match.(reg, filter(x -> occursin(reg, x), mdp_file)), 1))...)
    Name = Vector{String}()
    Natm = Vector{Int32}()

    locaters = findall(x -> occursin(r"\s*moltype.+\(", x), mdp_file)
    resID = zeros(Int32, length(Nmol), maxNatm)
    Catm = zeros(Int32, length(Nmol), maxNatm)
    Ratm = zeros(Float64, length(Nmol), maxNatm)
    Satm = zeros(Float64, length(Nmol), maxNatm)
    Eatm = zeros(Float64, length(Nmol), maxNatm)
    Qatm = zeros(Float64, length(Nmol), maxNatm)
    Tatm = Array{String}(undef, length(Nmol), maxNatm)
    for locater in locaters
        Imol = parse(Int32, match(r"\s*moltype.+\((\d+)\)", mdp_file[locater]).captures[1]) + 1
        push!(Name, match(r"\"(.*)\"", mdp_file[locater + 1]).captures[1])
        push!(Natm, parse(Int32, match(r"\((.*)\)", mdp_file[locater + 3]).captures[1]))
        locater += 3
        for i = 1:Natm[Imol]
            # 提取每个原子的类型信息
            txt = mdp_file[locater + i]
            idx = parse(Int32, match(r".*type=\s*(\d+).*", txt).captures[1])
            resID[Imol, i] = parse(Int32, match(r".*resind=\s*(\d+).*", txt).captures[1]) + 1
            Catm[Imol, i] = idx
            Ratm[Imol, i] = Rad[idx + 1]
            Satm[Imol, i] = sigma[idx + 1]
            Eatm[Imol, i] = epsilon[idx + 1]
            Qatm[Imol, i] = parse(Float64, match(r".*q=\s*([^,]+)\,.*", txt).captures[1])
        end
        # 提取每个原子的原子名信息, 在atom (原子数)后面
        locater += Natm[Imol] + 1
        for i = 1:Natm[Imol]
            Tatm[Imol, i] = match(r"name=\"(.*)\"", mdp_file[locater + i]).captures[1]
        end
    end

    # 提取每个残基信息
    locaters = findall(x -> occursin(r"\s*residue \(", x), mdp_file)
    maxNres = max(parse.(Int32, getindex.(match.(r"residue \((\d+)\)", mdp_file[locaters]), 1))...)
    resName = Array{String}(undef, length(Nmol), maxNres)
    Imol = 0
    for locater in locaters
        Imol += 1
        Nres = parse(Int32, match(r"residue \((\d+)\)", mdp_file[locater]).captures[1])
        for i in 1:Nres
            nr = parse(Int32, match(r".*nr=(\d+).*", mdp_file[locater + i]).captures[1])
            name = match(r".*name=\"(.+)\".*", mdp_file[locater + i]).captures[1]
            resName[Imol, i] = @sprintf("%05d%s", nr, name)
        end
    end

    # 从键角信息判断H的连接原子并指定H的类型
    locaters = findall(x -> occursin(r"^ +Angle:", x), mdp_file)
    for (Imol, locater) in enumerate(locaters)
        nangles = parse(Int32, split(mdp_file[locater + 1])[2])
        if nangles > 0
            for line in mdp_file[locater + 3:locater + 2 + nangles]
                if !occursin("ANGLES", line)
                    break
                end
                paras = split(line)
                k = parse(Int32, paras[end])
                j = parse(Int32, paras[end - 1])
                i = parse(Int32, paras[end - 2])
                if uppercase(Tatm[Imol, i + 1])[1:1] == "H"
                    Tatm[Imol, i + 1] = @sprintf("H%s", Tatm[Imol, j + 1])
                end
                if uppercase(Tatm[Imol, k + 1])[1:1] == "H"
                    Tatm[Imol, k + 1] = @sprintf("H%s", Tatm[Imol, j + 1])
                end
            end
        end
    end

    # 输出到qrv文件
    Ntot = 0
    Nidx = 0
    for i = 1:Ntyp
        for n = 1:Nmol[i]
            for j in 1:Natm[i]
                Ntot += 1
                if Ntot in ndxPro || Ntot in ndxLig
                    Nidx += 1
                    if radType == 0
                        radi = Ratm[j]
                    elseif radType == 1
                        radi = getRadi(match(r"([a-zA-Z]+)\d*", Tatm[i, j]).captures[1])
                    else
                        error_hint("Error: radType should only be 0 or 1. Check settings.\n")
                        readline()
                        exit(1)
                    end
                    write(qrv, @sprintf("%6d %9.5f %9.6f %6d %9.6f %9.6f %6d %s %s %-6s  ", 
                    Nidx, Qatm[i, j], radi, Catm[i, j], Satm[i, j], 
                    Eatm[i, j], Ntot, "\"$(Name[i])\"-1.$j", resName[i, resID[i, j]], Tatm[i, j]))
                    if Ntot in ndxPro
                        write(qrv, "Pro\n")
                    elseif Ntot in ndxLig
                        write(qrv, "Lig\n")
                    end
                end
            end
        end
    end

    close(qrv)
end

function getRadi(attype)# mBondi from AMBER20/parmed/tools/changeradii.py
    radBondi = Dict{AbstractString, Float64}(["C" => 1.7, "H" => 1.2, "N" => 1.55, "HC"=> 1.3, "O" => 1.5, "HN"=> 1.3, 
    "F" => 1.5, "HP"=> 1.3, "SI"=> 2.1, "HO"=> 0.8, "P" => 1.85, "HS"=> 0.8, "S" => 1.8, "CL"=> 1.7, "BR"=> 1.85, "I" => 1.98])
    attype = uppercase(attype)
    if length(attype) >= 2
        return get(radBondi, attype[1:2], radBondi[attype[1:1]])
    else
        return get(radBondi, attype, 1.5)
    end
end
